The interplay between screening properties and colloid anisotropy: towards a reliable 

pair potential for disc-like charged particles 
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The electrostatic potential of a highly charged disc (clay platelet) in an electrolyte is investigated 
in detail. The corresponding non-linear Poisson-Boltzmann (PB) equation is solved numerically, 
and we show that the far-field behaviour (relevant for colloidal interactions in dilute suspensions) 
is exactly that obtained within linearized PB theory, with the surface boundary condition of a 
uniform potential. The latter linear problem is solved by a new semi-analytical procedure and both 
the potential amplitude (quantified by an effective charge) and potential anisotropy coincide closely 
within PB and linearized PB, provided the disc bare charge is high enough. This anisotropy remains 
at all scales; it is encoded in a function that may vary over several orders of magnitude depending 
on the azimuthal angle under which the disc is seen. The results allow to construct a pair potential 
for discs interaction, that is strongly orientation dependent. 
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I. INTRODUCTION 

Clays, in the generic form of charged platelets, enjoy 
widespread use in applications ranging from drilling, rhe- 
ology modification (for paints, cosmetics, cleansers...), 
catalysis etc. As a significant component of soils, clays 
are also of importance for crop production. The difficulty 
of synthesizing clays with well controlled properties (size, 
composition, charge...) has long hindered their funda- 
mental study. The situation has considerably changed 
in the last ten years, with the increasing availability of 
customized synthetic clays, among which Laponite is a 
prominent example. Yet, our understanding of such sys- 
tems is rudimentary (see e.g. 0, S 0, 0, UIE HI an d 
references therein). 

A reasonable model for Laponite platelets is that of 
uniformly charged and infinitely thin discs 0. In this 
paper, the focus will be on electrostatic interactions be- 
tween identical charged discs, a crucial ingredient for un- 
derstanding the phase behaviour and stability of clays in 
suspensions. The high anisotropy of these objects makes 
analytical progress difficult. In addition, these discs are 
typically highly charged, and the electrostatic coupling 
with their electrolytic environment (microscopic charged 
species) needs to be described by non-linear theories: 
the plain linear Debye-Hiickel approach should fail. We 
will work here in the common framework of non-linear 
Poisson-Boltzmann (PB) theory, where the (dimension- 
less) electrostatic potential outside the charged macro- 



ions obeys an equation of the form V 



' sinh i 



assuming for simplicity monovalent microions only, the 
density of which governs the screening length k _1 (the 
Debye length). 

In a solution, the typical distance between macroions 
is often larger than the Debye length (this condition re- 
quires a minimal but nevertheless small amount of salt). 
At these "large" scales, the potential created by a given 
disc is small enough -compared to thermal agitation- to 
allow for the linearization of PB equation: V 2 </> ~ K 2 (f>. 
Accordingly, the potentials within non-linear PB on the 



one hand, and linearized PB theory with a suitably cho- 
sen boundary condition on the other hand, coincide at 
large enough distance from the colloids, be they of dis- 
cotic or other shapes. An analytical treatment within 
linearized PB (LPB) is of course considerably simpler 
than within PB, but the above remark may be of little 
practical help if one is not able to derive the relevant 
boundary condition on the colloid (effective potential), 
such that the corresponding LPB solution reproduces the 
PB one in the region of low enough potential. Close to 
the colloids, non-linear effects prevail (LPB and PB so- 
lutions strongly differ), and broadly speaking, microions 
-essentially counterions- suffer there a high electrostatic 
coupling and may be considered as "bound". They de- 
crease the bare charge of the colloid so that its elec- 
trostatic signature at large distance defines an effective 
charge which is usually smaller in absolute value than 
the bare one (close to the colloid, the effective potential 
is accordingly smaller than its non-linear counterpart). 

For a unique charged sphere in an electrolyte, PB and 
LPB theories give rise to the same far- field behaviour, of 
Yukawa type [exp(— nr)/r, where is r is the radial coor- 
dinate]; non- linear effects only affect the prefactor (from 
which the effective charge is defined) , preserving the func- 
tional form of the potential [ToLfTTl 121. The same remark 
equally applies to an infinite rod. The situation changes, 
however, for anisotropic objects such as discs or finite size 
rods [l3j |. where non- linear screening phenomena generi- 
cally affect the functional form of the potential and can- 
not be subsumed in an effective scalar quantity (effective 
potential or charge). In other words, whereas the symme- 
try of the effective colloid clearly remains spherical in the 
case of spheres, predicting the symmetry of the effective 
charge distribution and associated electrostatic potential 
for a highly charged disc is a non trivial question. 

It is the purpose of the present work to study how non- 
linear screening effects and anisotropy conspire to affect 
the far-field behaviour in the case of discs. It has been 
shown in that highly charged spheres and infinite 
rods may be considered as objects of constant effective 
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potential e ff (in the sense that 4>eS becomes independent 
of physico-chemical parameters, provided that na > 1 
where a is the colloid radius. The complementary results 
reported here indicate that the constant potential pic- 
ture goes in fact beyond this analysis, and give the cor- 
rect symmetry of the effective charge distribution onto 
the disc. Such a boundary condition (within LPB the- 
ory) produces the same electrostatic potential as a highly 
charged disc within PB. A physical argument allowing to 
anticipate this correspondence will be presented in sec- 
tion [H] Since the exact LPB solution for a disc held 
at constant potential in an electrolyte is not known, we 
will introduce in section IIIII a semi-analytical procedure 
to solve this problem. The characteristic features of the 
electrostatic potential relevant for clay discs will be ob- 
tained, and in order to assess the validity of the constant 
potential picture, the corresponding electrostatic poten- 
tial will be compared in section IIVI to the solutions of 
the full non-linear PB theory. The latter will be ob- 
tained through an iterative numerical procedure. From 
these results, a pair potential will be constructed for 
charged discs, that has the same status as the celebrated 
Derjaguin-Landau-Verwey-Overbeek expression relevant 
for spheres [H [H El , and which includes charge renor- 
malization. Concluding remarks will finally be presented 
in sections IVl and IVT1 



II. THE CONSTANT EFFECTIVE POTENTIAL 
PICTURE: WHY? 

A disc of radius a and uniform surface charge den- 
sity CTbarc = -^barc e/(7ra 2 ) is immersed in an infinite 
sea of electrolyte with bulk density n s . From the per- 
mittivity e of the solvent, the Bjerrum length is de- 
fined as £b = e 2 /(ekT), where kT is the thermal energy 
and e denotes the elementary charge. Within non-linear 
Poisson-Boltzmann theory, the dimensionless potential 
\4> = eip/(kT), ip being the original electrostatic poten- 
tial] created by the disc obeys the following equation 



V 2 = n 2 sinh(0) 



(1) 



where k is the inverse Debye length defined through k 2 = 
87r£sn s . For convenience, </> is chosen to vanish at infinity. 
Eq. holds outside the disc. 
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Considering a highly charged disc with furthermore 
Ka > 1, one may partition space into three regions, as 
sketched in Fig. Q In region A, the electrostatic coupling 
between the colloidal disc and the microions is most im- 
portant and one has (j) > 1. Outside A, in regions B and 
C, one has (j> < 1 and PB equation Q may be linearized: 
the corresponding Helmholtz-like LPB equation reads 



V 2 c 



K 2 (j). 



(2) 



In addition, in region B, cj) is of unidimcnsional character 
and well approximated by the potential created by an in- 
finite charged plane. A contrario, in region C, <j) regains 
its full 3D nature (2D here with the present azimuthal 
symmetry). The lateral extension of the "non-linear re- 
gion" A is given by k _1 while a measures the extension 
of region B. Since we assume na > 1, we have A C B and 
moving away from the disc, non-linear effects disappear 
before the finite size of the disc becomes relevant. In 
other words, B is the non vanishing intersection between 
the "linear" region and its one-dimensional counterpart 
where the solution of Eq. takes the form ^| 



iiD = 4arctanh(7e KZ ) 
~ 4je~ KZ in region B. 



(3) 
(4) 



In these expressions, z denotes the distance to the plane 
and, assuming without loss of generality a positive bare 
charge CTbare , 7 is the positive root of the quadratic equa- 
tion 



7 



7 + 1 = 0. 



(5) 



FIG. 1: Schematic side view of the charged platelet 



In colloidal dispersions, the relevant range for the inter- 
actions is that of far-field (except for dense systems) and 
the behaviour in the non-linear region (A) is of little in- 
terest. From Eq. 0}, it appears that the potential felt 
in the linear region B+C, when extrapolated to contact 
(z = 0), reads 4> c s = 47. As a consequence, solving LPB 
equation (j2Jl with this boundary condition should provide 
the same potential outside region A as the PB solution 
of equation 0J. 

The above remarks follow from the constraint a > k^ 1 
and apply irrespective of the value of the bare charge. 
In particular, CTbarc (hence e ff through 7) may be posi- 
tion dependent on the disc. However, we are interested 
here in highly charged discs for which the non-linear re- 
gion A exists (for low CTbare, region A disappears; PB 
and LPB solutions coincide at all distances and the is- 
sue of effective potentials becomes trivial: effective and 
bare charges are equal). From Eq. JSJ, it appears that 
< 7 < 1 and that 7 — > 1 when CTbare becomes large, so 
that 4> c f{ — > 4, and the field created is independent of the 
bare charge. More details concerning the phenomenon of 
effective charge saturation may be found in 

We conclude here that a highly charged disc should 
effectively behave as a constant potential object treated 
within a linear theory. The corresponding LPB prob- 
lem will be addressed in the following section but we 
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emphasize before that the argument developed here pro- 
vides 4> e g to leading order in na. On the basis of the 
behaviour of charged spheres for which the curvature cor- 
rection has been computed in |16| [leading to the result 
4> e & = 4 + 2/(1 + na) + 0(ko,)~ 2 }, we anticipate that </> e ff 
may exceed the threshold 4. A similar behaviour is ob- 
served for cylinders of infinite length |l6j . 



III. SEMI-ANALYTICAL SOLUTION OF THE 
DIRICHLET LINEARIZED PB PROBLEM 

A. Methodology 

If the solution of LPB equation © with Dirichlet 
boundary condition (j> = (f)Q was known analytically, the 
effective charge of highly charged discs would follow 
immediately, enforcing <f>o — 4 on the disc surface. Unfor- 
tunately, such a solution only exists in vacuum (i.e. when 
k = 18]). To our knowledge, the only solution known 
at finite k is that associated to Neumann boundary condi- 
tions (uniform surface charge) 0] . To solve the Dirichlet 
problem, we have therefore developed a semi-analytical 
procedure where the problem at hand is recast into a 
Fredholm integral equation (see below and appendix lA^I . 

The general solution of Eq. J2J) may be written assum- 
ing both cylindrical symmetry around an axis (Oz), and 
reflection symmetry z <-> — z: 

/>oo 

ct>(p,z)= A(k)J (kp)e- Vk ' A+K ^dk. (6) 
Jo 

In this relation, (p, z) denote the cylindrical coordinates 
and Jo is the Bessel function of order 0. The difficulty 
in the present situation is that the boundary conditions 
imply that the weight function A(k) obeys the mixed 
system 

/>oo 

/ A(k) J (kp)dk = 0o for p < a (7) 
Jo 

/>oo 

/ y/W^- A(h) J (kp) dk = for p > a (8) 
Jo 

where the second equation follows from the vanishing 
of the normal electric field d z 4> on the symmetry plane 
z = 0. Starting from Eqs. JJJ) and (jHJ), the problem 
is rephrased in terms of an integral equation, solved nu- 
merically, from which the function A(k) is computed (see 
Appendix^}. The potential then follows from ©. 

B. Properties of the solutions 

A typical solution is shown in Fig. [5] In the remainder, 
the variable 9 <E [0, 7r/2] denotes the angle between a 
given direction and the normal to the disc (9 — n/2 in 
the symmetry plane z = and 9 — along the normal to 
the disc, i.e. when p = 0). It may be observed that the 




FIG. 2: Solution of LPB equation with Dirichlet bound- 
ary condition of a uniform potential <f> — <f)o on the disc sur- 
face. The quantity rcj) is shown on a linear-log scale to empha- 
size the far-field behaviour. Here, na = 2 and r = sj p 2 + z 1 
denotes the distance to the disc center. The continuous curve 
shows the potential in the direction 8 = n/2 (as a function of 
r/a — p/a), whereas the dashed line shows the behaviour as 
a function of r/a = z/a along normal axis p = (9 — 0). The 
inset shows <f>/4>o on a linear scale, again in the two perpen- 
dicular directions — and 6 — n/2. 



potential is anisotropic at all distances, a generic feature 
of screened electrostatics 0,^3 : the behaviours for 9 = 
and 9 = ir/2 strongly differ, at all scales. The anisotropy 
of the potential at large scales is encoded in a function 
f(Ka,9), such that expression © may be written, for 
nr > 1, 

-nr ( e ~ KT \ 

<P(r,6) ~ Zl B f( K a,e) +0 — , (9) 

r \ r / 

where r = (p 2 + z 2 ) 1 / 2 again denotes the distance to 
the disc center. In Eq. @, the total charge Z of the 
platelet appears. Z is the integral over the disc surface 
of the surface charge density a(p) [Ze — f <r(s)d 2 s\. 
This density turns out to be related to the anisotropy 
function through 

f{na,0)= [ ^ exp(-/cf • s) d 2 s, (10) 

./disc % e 

f being a unit vector pointing in the 9 direction. As 
expected, without electrolyte, k vanishes so that / = 1 
and the potential in © takes the familiar form of an 
isotropic Yukawa expression. 

The anisotropy function / and the total charge Z, are 
the key quantities governing far-field behaviour. Note 
that Z is not known a priori since only the surface po- 
tential is imposed. It may be shown that / is related to 
the weight function A(k) appearing in JBJ through 

Z e „. „, Aii Ksin#) 

_/<„,*) = -_L_^. (ii) 
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For small arguments, one has A(x) oc x so that /(««, 0) = 
1, which also means that the total charge is directly ac- 
cessible through the behaviour along the 9 = axis : 



0(r,O) 



(12) 



Note also that f(na, 0) = 1 directly follows from i|10|) 
since f • s = when 9 = 0. 




Ka=l 
Ka=2 
Ka=4 - 



JT/4 

e 



Jt/2 



FIG. 3: Anisotropy function for kcl=1, 2 and 4 as a function 
of azimuthal angle. Fhe inset shows the results on a linear-log 
scale, for na = 5 and Ka = 10. 

Figurc[3]shows /(«a, 9) for different salinity conditions. 
This function increases with 9 and may take large values 
when Ka exceeds a few units (see the inset where the 
y-axis is shown in log scale). On the other hand, for 
Ka < 1, / remains close to unity in all directions. The 
potential is strongest in the disc plane (8 = 7r/2), and 
increasing screening (k), one also strongly increases the 
anisotropy of the electrostatic potential. For a reasonable 
value Ka = 10, f(9) varies by almost 3 orders of magni- 
tude (a factor 930). In figure^] the complementary infor- 
mation concerning the total charge Z is displayed. This 
quantity will be further discussed in section IlII CI Note 
that /(/ca, 9) does not depend on <po, since it probes the 
repartition of surface charge distribution, not its overall 
magnitude. On the other hand, Z scales linearly with (f>Q. 



C. An approximate expression for the anisotropy 
function and charge 

It is instructive and useful for practical purposes to 
have an approximate analytical expression for the func- 
tion f(Ka,6). From Eq. I|10|) . that may rewritten 



f{Ka,9) = — J I {Kp sin 9) a (p) pdp, 



(13) 



this amounts to look for an approximate expression for 
the surface charge density a. To this aim, we recall [2(j 
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8 10 



FIG. 4: Charge 4Z(.B/(a<j>o) as a function of na (continuous 
curve). The dashed line shows the result obtained within the 
simplified two parameter model [Eqs. 1)15^ and p7^ - ^18| . see 
text]. The dotted line corresponds to Eq. 1201 . 



that for an ideal conducting disc when k = 0, a diverges 
in the vicinity of the edge (p — > a) since 



e^o 1 

2TT 2 a£ B y/1 - (p/af 



(14) 



Recall that (j)Q denotes the dimensionless electrostatic po- 
tential = ip/(ekT). The singularity of the electric field 
near sharp edges, which is the reason for the efficiency of 
lightning conductors, also pertains in presence of an elec- 
trolyte. We indeed show in appendix iBl that when k ^ 0, 
the Dirichlct solution to Eq. (J2J exhibits a similar diver- 
gence as that present in Eq. (|14f> . namely a oc (a— p) -1 / 2 . 




FIG. 5: Electrostatic potential following from the two pa- 
rameters approximation (circles) compared to the exact result 
(continuous curve, obtained following the method detailed in 
section ITTrAl and Appendix|X). The main graph corresponds 
to 9 = and the inset to 9 = n/2. 

A very simple two parameter ansatz fulfilling this di- 
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FIG. 6: Anisotropy function resulting from the two parame- 
ters approximation 115H where oo and <ti are determined from 
l|17fl - l|18|L compared to the exact result. The two above pa- 
rameters are plotted in Fig. |7| Here, na — 5 whereas na — 1 
in the inset. 




Ka 



FIG. 7: Dimensionless charges a' 01 = 2na(.B 0"o,i/(e0o) fol- 
lowing from p5p supplemented with Eqs. I17p - Q18| . as a func- 
tion of salinity conditions. 



vergence requirement is 



<r(p) 



2 VI- (p/af 



(15) 



From this expression, the anisotropy function may be 
computed and takes the form 



Co 2ii(/casin( 



<j\ sinh(Kasin( 



f( K a,6) = 

u o T w 1 /^ u sin [/ c o "T i ft a sin 1/ 

(16) 

In Eqs. i|T3|) and JTBJ, Jo and /i denote modified Bessel 
functions of the first kind, of order and 1. Expressions 
(fT5*|) and (flf))) are not exact and there are several ways 
to choose the two parameters Co and <7i, that will be 
determined by two constraints. The simplest possibility 



is to enforce 0(0, 0) = 0(a, 0) 
the choice (hereafter adopted) 



0(0,0) = O 
(0(p,O)) = 0o 



but it turned that 

(17) 
(18) 



gave better results (the angular brackets denote average 
over the disc surface). In the limit k — > 0, expressions 
(|15|l and (|16fl become exact (with ctq =0), and we expect 
that the comparison with exact results at finite k will be 
all the better as na is low. It is indeed the case, but when 
Ka = 5, the approximation is still reasonable (see Figs. 
[S] and EJ) . By comparison with the exact solution Fig. 
shows that the anisotropy of the potential is correctly 
captured. The corresponding values of partial surface 
charges <7o and a\ are shown in Fig. where it may be 
observed that <7o vanishes at low salt, as expected. The 
associated total charge Z is given by 



Z 



(oo + en ) 



O0o 

2£r 



(a'a+a'J, (19) 



where <t[ and a' 2 are the quantities plotted in Fig. 
The charge Z is displayed in Fig. 0] and compares favor- 
ably with its exact counterpart. In the above expression, 
however, oo and a\ are functions of na with unknown 
analytical expression. As it is desirable to have an an- 
alytical formula, we propose the following argument: in 
the limit of large na, the disc essentially behaves as a an 
infinite plane, from which we deduce ZIb/cl ~ 0o«a/2. 
To estimate the next order correction [constant term C 
in the expansion Z£b/(ckPo) = {^Ka/ + C)/4] we may 
take the limit Ka — where the solution is given by {TIJ, 
which imposes ZIb/cl = 0o/vr. We therefore obtain 



z ^b_ _ 00 
a ~ 4 



2 Ka + — 

7T 



(20) 



Anticipating that the relevant values of 0o are close to 
4 (see sections |n] and II V B , we have factorized the ratio 
0o/4 in the previous relation. The quality of approxi- 
mation (|20|l is assessed in Figure 0J which shows a good 
agreement. On the other hand, extracting the correction 
factor C from the large Ka behaviour of the exact Z dis- 
played in Fig. 0] gives C ~ 1.88, to be compared with 
C = 4/tt ~ 1.27 in Eq. {ZjJ. 



IV. NUMERICAL RESOLUTION OF THE 
NON-LINEAR POISSON-BOLTZMANN THEORY 



In section IIIII we have obtained the solution of lin- 
earized PB theory with uniform potential boundary con- 
dition on the disc. From the discussion developed in sec- 
tion HH we expect the properties described (with 0o = 4) 
to characterize also the far-field created by a highly 
charged disc, treated within non-linear PB theory. In the 
following, we critically test this scenario. We first present 
the numerical procedure used to solve the non-linear PB 
problem. 



6 



A. Green's function formalism and numerical 
method 

Introducing explictly the charge density q d (r) borne by 
the disc, Eq. Q is rewritten 

V 2 = K 2 sinh(^)+4^ B ^^. (21) 

e 

In the subsequent analysis, we will consider the case of a 
uniformly charged disc for which one has, in cylindrical 
coordinates 



(|27|l ]. In addition, it seems that the optimal choice for 
the numerical screening parameter kq is kq ~ n. We have 
checked that the solutions found were independent of n 
(as they should) by changing this parameter in the range 
[k/5; 5k]. The previous procedure bears similarities with 
the one used in Ref. |2lj . where a confined geometry 
was considered (whereas the situation is that of infinite 
dilution here). 



B. Results 



q d (r) = a harc S(z)Q(a - p), 



(22) 



where O is the Heaviside step function and 5 the Dirac 
distribution. However, it is important to emphasize that 
the results that will be derived are more general, and 
hold irrespective of the precise PB boundary condition 
on the disc, provided the bare disc charge is high enough 
(phenomenon of effective potential saturation). 

In view of a numerical resolution, it is convenient to 
rewrite (|21|l in the form 



(V 2 



K 2 sinh(0) + 47r^> gd (r)/e-K 2 0, (23) 



where kq is an arbitrary quantity that will be optimized 
in order to speed up the resulting procedure (see below) . 
Introducing the Green's function 



„-Ko|r— r'| 

47r|r - r'l 



solution of 

(V 2 -K 2 )C?(r,r') 
we may recast into 



6(1 



(r)=|e(r,r') 



k 2 sinh[</>(r')] + 4mli 



fir' 

The contribution arising from q d may be computed ana- 
lytically: 



g(r,r')q d (r')dr' 



ZL 



J 1 K)J (M eXP( 1 Z '^ T ^ ) rffc(27) 



and Eq. I|2t)|) is solved iteratively. Starting from the ini- 
tial guess 0o = 0, the right hand side of l|2tj[l . denoted 
(j>Q nt is computed. This provides a new input potential 
4>i = a<f>^ ut + (1 — a)4>o, which is itself inserted in the 
rhs of 112611 to produce <f>° ut etc. The mixing parame- 
ter a is chosen in the range [10~ 2 ; 1CP 1 ]. Convergence 
4> n ~ 0° ut is generally achieved for typically 50 to 200 
iterations. The procedure may be accelerated starting 
not from 0o = but from the solution of LPB theory 
[known in the present Neumann case, and given by Eq. 



From the method sketched in section llV Al the numer- 
ical solution of the non-linear PB equation Q may be 
obtained for arbitrary bare charges and salt content. 




(24) 



FIG. 8: Dashed line: electrostatic potential solution of PB 
equation Jy for a highly charged disc with uniform bare sur- 
(25) face charge (Neumann-like boundary condition, Z^ arc £B /a = 
60). Here sa = 5. The continuous line shows the solution of 
LPB theory with Dirichlet boundary condition 0o = 4 (see 
section HTH . The behaviour is shown along the two perpen- 
dicular directions — and 9 = tt/2. The main graph shows 
rij) versus r/a = (p 2 + z 2 ) 1 ^ 2 / 'a on a linear-log scale. The inset 
shows the previous potentials vs r on a linear scale. 



To test the constant potential picture put forward in 
section [H] (which dwells on the fact that na is "large 
enough") we show in Fig. [S]the PB potential correspond- 
ing to a "large" bare charge. It appears that the PB and 
LPB potential are in excellent agreement except in the 
immediate vicinity of the disc, so that the constant effec- 
tive potential prescription seems accurate. 

From the PB potentials, we may also extract the ef- 
fective charge Z c g and anisotropy function f(na, 9), that 
convey a more complete information than a plot like that 
of Fig. |S| Z c s follows from the far-field behaviour along 
the = axis [see Eq. (TJJ}]: 



(r,0) ~ Zeff^i 



Once Z e g is known, / is computed from Eq. i jffi l. 



(28) 
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FIG. 9: Effective charge as a function of the bare one, for 
na — 5. The circles correspond to the PB solution, and the 
dotted line has slope one to show the regime of weak coupling 
where Z e g = Z bmQ . The analytical expressions for a plane, a 
cylinder and a sphere have also been plotted (see labels). 
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FIG. 10: Saturated effective charge Z^g as a function of 
salt (circles). The LPB charge Z is also shown for <j>o = 
4 (continuous line). The dashed line displays the empirical 
expression 1291 . 



For fixed Ka, Z c s is a function of the bare charge, 
see Fig. [5J where the corresponding analytical predic- 
tions for planes, spheres and cylinders (of infinite length) 
have also been reported 0. Not surprisingly, the be- 
haviour is intermediate between that of infinite planes 
and spheres, and somehow resembles the results valid 
for charged cylinders. In addition, Z e g reaches a satura- 
tion plateau when Zb arc becomes large Q (see Fig. [§J. 
This asymptotic plateau defines the effective charge at 
saturation Z^g, which is shown in Fig. El This quan- 
tity is an increasing function of salt content (except for 
Ka < 0.5, see below), since an increase in salt density 
enhances the macroion/microion screening, which dimin- 
ishes the amount of "counterion condensation" and con- 
sequently increases the effective charge 0] . Figure El 



shows that the constant potential prescription of section 
[njwith cj)Q = 4 provides a satisfying description of highly 
charged platelets, as far as the (effective) charge is con- 
cerned. It also appears that the best linear interpolation 
reads 



Zl$£ B /a~ 2/ta + 2.9 



(29) 



which is rather close to approximation (|20|l with the 
choice 0o = 4. It may be observed in Fig. El that for 
low Ka, the effective charge increases when Ka decreases. 
This aspect will be discussed in section IV Bl 

Computation of anisotropy functions confirms the rel- 
evance of the constant effective potential picture (see the 
comparison proposed in Fig. Illfl . As shown in the inset, 
the very high values / ~ 900 predicted from the analysis 
of section lllll are indeed found within non-linear PB. Note 
that the agreement reported in Fig. 1111 is only expected 
at high bare charges. For low bare charges, f(Ka,8) de- 
pends -at variance with its large bare charge counterpart 
-on the details on the boundary conditions chosen on the 
disc to solve PB. In the present situation (uniform sur- 
face charge), / may be computed analytically with the 
result H 



f(Ka,8) = 2 



I\(Ka sin#) 



Ka sine 



(30) 



This functional form is observed from our numerical data, 
for Z bm . c £ B /a < 1 (see Fig. IT2")) . It turns out to differ 
much from that reported in section IIIII (shown with a 
dashed line in Fig. 11211 . We may also observe in Fig. El 
that the "Neumann" expression (|30|l is lower than the 
Dirichlct one. The reason is the following : / is sensitive 
to the charges lying near the edge of the disc [see e.g. 
Eqs. (|10|) and (|13fl ]. With Dirichlet boundary condition, 
the induced surface charge diverges (see appendix 0), at 
variance with the situation of a uniform surface charge, 
which therefore exhibits a less anisotropic potential. The 
Dirichlet and Neumann expressions respectively provide 
upper and lower bounds for the anisotropy function. 

In spite of the good agreement shown in Fig. [SJ a 
slight difference may be observed between PB and LPB 
results. It may be concluded that a value (/>o slightly 
above 4 may give a better agreement between non-linear 
and linear profiles. As mentioned at the end of section 
ITT1 finite Ka corrections increase the value of the effective 
potential above 4 for spherical and rod-like macroions. 
Figure El an d El show that a similar effect exists for 
discs : when the 0o of LPB approach is considered as an 
adjustable parameter, the agreement between PB and 
LPB (again for highly charged discs) becomes excellent 
even at relatively low values of Ka, such as Ka = 0.5 
or even Ka — 0.1. The salt dependence of the above 
optimal potential, denoted 4>r] pt , is shown in Fig. El 
One observes that 0Q Pt is close to 4 for Ka > 5, but 
may take significantly different values at lower Ka. This 
leads to reconsider the plot of Fig. El smce the (relative) 
discrepancy PB /LPB may arise from discarding finite Ka 
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B=k/2 



p/a 



FIG. 11: Anisotropy function f(na,0) as a function of az- 
imuthal angle, for so = 5. The crosses correspond to the PB 
result with Z^ avc £b /a — 60 (i.e. in the saturation regime) 
while the continuous curve shows the constant potential LPB 
result obtained in section ITTT1 Inset: same for na — 10. 



L,8 r 



Ka=2 




FIG. 12: Same as Fig. 1111 for a low bare charge Zbare (-h/cl 
10" " 
0.5 

continuous line. For the sake of comparison, we also plot 
with a dashed line the constant potential LPB result already 
displayed in Fig. 1111 which is relevant at high bare charges. 



2 . The circles, squares and triangles correspond to ko = 
2 and 5. The prediction of Eq. 1301 is shown by the 



FIG. 13: Plot of the PB potential (continuous line) versus 
distance from disc center in the 9 = n/2 direction (disc plane) 
Also displayed are the LPB results for (j> = 
cles) and for </>o — 4 (dotted curve). Here, na 
-^bare^-B /a = 15 corresponding to the saturation plateau (the 
precise value of ^barc is therefore irrelevant) . The inset shows 
the same quantities along the disc normal. 



° pt = 5.7 (cir- 
0.5 and 




FIG. 14: Same as Fig. Uglfor Ka = 0.1, with now <j) = ^ pt 
9.0. 



uation of two discs (radii a\ and 02) in an electrolyte 



effects (i.e. enforcing 0o = 4). Fig. ^| compares PB 
saturated effective charge to its LPB counterpart, with 
4>0 = 4>o Pt ■ Both quantities now agree very well. This 
latter comparison is a severe and successful test for the 
relevance of the constant potential prescription. 



V. DISCUSSION 

A. Pair potential 

From the properties of the one body electrostatic po- 
tential (f> discussed previously, one may obtain the large 
distance behaviour for the pair potential f 12 in the sit- 



U12 = Z cSi1 Z cSi2 J(k a 2 ,0 2 ) ■ (31) 

r 

Here 0i is the angle between the normal to disc i and 
the center-to-center direction r 12 (with |r 12 | = r). The 
validity of such an expression at intermediate or short 
distances (i.e. nr of order 1) is unclear, since polarization 
effects of disc i on disc j should at least perturb the 
symmetry of the effective charge distribution, and hence 
alter the one body expression for / plotted in Fig. ^2 
We also note that the sub-leading terms in Eq. (|3T|) 
that become more important as nr decreases, involve a 
more complex dependence on relative orientations (with 
all Euler angles becoming relevant, contrary to the far- 
field case where only Q\ and 82 matter). 
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(> 

4 - - 



2 4 6 8 10 

FIG. 15: Circles: Values of the optimal surface potential 
</>o Pt to be imposed within LPB theory to produce the same 
far-field behaviour as a highly charged disc treated within PB 
(saturation limit). The dashed line indicates the high salt 
limiting behaviour (j>o = 4. 




FIG. 16: Same as Fig. EU but with <j) = <#J pt (the latter 
quantity being plotted in Fig. I15H instead of <j>o = 4. For 
completeness, the results obtained with cf>o — 4 are shown by 
the dotted line. 

We may conclude here that at fixed center-to-ccntcr 
distance r, the favored configuration is that where the 
discs are parallel and perpendicular to their ccntcr-to- 
center vector ri2. The T-shape configuration is interme- 
diate and the most repulsive one corresponds to coplanar 
discs (parallel to ri2, as coins lying on a table). How- 
ever, the situation changes if one fixes the closest distance 
D between the two discs. Comparing the configuration 
9\ = 02 = 7r/2 where T> — r — a\ — 02 to that with 
61 = O2 = for which T> — r requires to compare Q 
defined as 

Q = f (*Oi, ~) f (Koa, |) e-^+ a ^ (32) 

with 1. From approximate expression i|30|) . it appears 
that Q is always smaller than 1. Instead of l|30|) . a more 
reliable expression for the anisotropy is provided by (ll(j|) , 



which leads to the same conclusion. We therefore recover 
the intuitive result that the less repulsive configuration 
at fixed V is for Q\ =62 = tt/2 (two coins on a table). 

B. Behaviour at low na 

In the present study, we have focused on the regime 
na > 1 since according to the argument of section [HI 
it corresponds to the situation where the effective po- 
tential may be predicted analytically. It appears that 
the Debye length acts as a local probe to reveal the 
anisotropy of the macroion under study. Hence, in the 
limit of small na where this probe cannot resolve the 
disc dimension, we found that the anisotropy disappears: 
f(na = 0,8) = 1, V0. We may then speculate that at 
small na, the precise form of the macroion becomes irrel- 
evant so that we should recover the same results as for 
spheres [23| . From the analysis of Ramanathan [24[ , we 
may consequently expect in the saturation regime : 

ZS?— K - <1 -21n(«a) + 21n[-ln(«o)]+41n2 (33) 
a 

Such an expression diverges for na — > 0, indicating that 
potential (or charge) renormalization ultimately becomes 
irrelevant. However, with the lowest value of na in- 
vestigated in this work (na = 0.1), we have measured 
~ 6.1 a /£ B (see Fig. EJ), whereas Eq. Q gives 
approximately a value 6.7. Whether this agreement is 
incidental or not is unclear. We also note that if the 
discs behave as spheres for very low na, their saturated 
effective potential should coincide with Z^g Ib /a. This 
is not the case for na = 0.1 where we measured </>Q Pt ~ 9. 
This may indicate what a salinity condition na = 1/10 is 
not low enough to enter the "spherical" regime. 

C. Validity of the PB approach 

We now discuss the validity of the Poisson-Boltzmann 
theory underlying the present analysis. Such an ap- 
proach neglects microionic correlations (be they of elec- 
trostatic or other origin, such as excluded volume) 
while macroion/microion electrostatic correlations are 
correctly incorporated. In the vicinity of the charged 
discs where the counterion density may become large, the 
neglected correlations are most important, and may in- 
validate PB theory if the disc bare charge is too large (say 
■Zbarc > ^barc)- Since we have considered here the situa- 
tion of high Zbare to explore the PB saturation plateau, 
we need to justify the relevance of such a plateau. In 
other words, this amounts to elucidating the circum- 
stances under which Z sat < Z£are> since for ^>are > Zaat, 
one has, roughly speaking, Z eS ~ Z s&t . 

For a salt-free system, Netz has considered the valid- 
ity of PB theory in planar, cylindrical and spherical ge- 
ometries [25] . Since no general result exists, we present 
here a a simple argument concerning discs, which goes as 



10 



follows (sec 12] for the spherical case). Microionic corre- 
lations may be accounted for by the coupling parameter 
r = Is/d where d is a typical distance between microions 
in the double-layer. This distance is bounded from below 
by that, denoted d* , where all Z monovalent counterions 
are artificially condensed onto the disc (as would happen 
in the low temperature limit). It may be estimated writ- 
ing that the typical surface per microion (d*) 2 on the disc 
is the mean value ira 2 /Z. Hence, 



r ~ yJPgObnc ^barc^/(7ra 2 ). (34) 

For microions with valency z, one would obtain 



: 3 Zbarc4/(™ 2 



(35) 



For the sake of the argument, the factor it could be omit- 
ted. The important point here is that Zbare ^s/a may be 
large, which corresponds to the saturation regime of PB 
theory, with still T < 1, which justifies the mean- field 
assumption underlying PB. With typical Laponite pa- 
rameters and monovalent microions, we have T ~ 0.7 
for Zbare — 700, a reasonable value for the charge. In 
addition, Zb are ^s/a — 33 which is well beyond the linear 
regime where effective and bare parameters coincide (see 
Fig. El or Fig. 1171 corresponding to a lower salt concen- 
tration for which Zbare^s/a — 33 lies in the saturated 
region). We finally note that the intersection between 
the PB saturation regime and the consistency condition 
r < 1 is all the larger as a/is is big 26]. But of course, 
when Zb aro strictly diverges, T exceeds a few units and 
PB breaks down. Finally, to be specific, would be 

defined as the value of Zb&re such that the coupling pa- 
rameter r is of order unity. The case of Laponite seems 
somehow borderline, since .Zbare does not differ much 
from Z^™ T e (on the order of a few thousands) . 




D. The case of asymmetric electrolytes 

Bearing in mind the classical rule that increasing the 
valency of microions decreases the range of validity of 
PB theory [see Eq. j35jl]. we may extend previous re- 
sults to 1:2 and 2:1 electrolytes. The key ingredient in 
the constant effective potential prescription is indeed the 
analytical solution of the planar (ID) PB equation. The 
latter problem has been solved by Gouy almost a cen- 
tury ago |27| . and it turns out that the counterpart of 
the monovalent result 4>q = 4 reads <fio — 6 for 2:1 elec- 
trolytes (i.e monovalent counterions and divalent coions, 
or more precisely when the ratio of coion to counterion 
valency equals 2). In the reverse 1:2 situation, we have 



6(2-^3) ~ 1.608. 



(36) 



The 1:2 effective potential is smaller that the 2:1 po- 
tential since screening by monovalent counterions is less 
efficient than with divalent ones (hence a higher effective 
potential, and a higher effective charge 28] ). 

By simply plugging the above expressions for <j)Q into 
the expressions derived in the previous sections for sym- 
metric electrolytes, one may describe 1:2 and 2:1 situa- 
tions as well. We finally note that the electrolyte asym- 
metry does not affect the anisotropy function f(na,6): 
LPB equation takes the same form (modulo a change in 
the numerical value of k) and only (j)Q is affected. 



E. Comparison with existing results 

In Ref we have addressed a similar issue as in the 
present paper. However, neither the LPB at constant 
potential nor the PB theory were solved. The anisotropy 
function has been estimated there from the Neumann 
LPB result with a uniform surface charge. This leads 
to expression (|30|l , which has been shown in Figure IT51 
be an underestimation of the Dirichlet result [and the 
agreement between PB and Eq. (|30|l is simply due to the 
low charge in Fig. ^1 see Fig. 111) . 

Following similar lines, the saturated effective charge of 
discs have been estimated in £| from the Neumann LPB 
solution. When the resulting dimensionless potential is 
equated to 4 on the disc center, we obtain 4] 



2 na 



£b 1 — exp(— na) 



(37) 



bare b 



which is very close to 2na as soon as na > 2. Such 
an expression only captures the leading order behaviour 
(the planar limit), but misses the offset correction (2.9) 
as appears in Eq. 



VI. CONCLUSION 



FIG. 17: Same as Fig. |*j|with na = 1. 



We have presented in this paper a detailed compari- 
son between the electrostatic potentials obtained within 
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Poisson-Boltzmann (PB) and Linearized PB approxima- 
tions, for a charged disc in an electrolyte. We have pro- 
posed a new and efficient semi-analytical method to solve 
the LPB problem at constant surface potential 4>o- The 
procedure used is not restricted to the specific problem 
considered here, and allows to solve more general situa- 
tions of the form given by Eq. I|A3(1 . On the other hand, 
the PB problem has been solved numerically following 
similar lines as in Ref [2]]]. We have shown that the far- 
field potential created by a highly charged disc within PB 
is remarkably close to its LPB counterpart with a suit- 
ably chosen value of cj>o , which therefore defines the effec- 
tive potential of the charged discs. As expected from the 
argument put forward in section^ the latter quantity is 
close to 4 (meaning that the effective surface potential is 
close to AkT/e) whenever na is larger than a few units, 
say kcl > 3. These results extend the conclusion of Ref. 
|14|. The argument of section [H] should in fact apply for 
any charged macro-ion for which the curvature is smaller 
than the inverse Debye length n of the surrounding elec- 
trolyte. Note also that in the limit of high bare charge, 
the details of the bare charge distribution onto the discs 
are irrelevant. In this respect, the results obtained here 
within PB with uniform surface charge are generic and 
would resist to charge modulation. 

The scenario emerging is that due to non-linear screen- 
ing phenomena, highly charged macroions may be consid- 
ered as effective objects with a uniform surface potential 
and can be treated within a linear theory (provided short 
distance features are irrelevant) which considerably sim- 
plifies the analysis and opens simulation routes. In addi- 
tion, this potential is constant provided there is enough 
salt in the solution, in the sense that it no longer de- 
pends on physico-chemical parameters. Such a viewpoint 
not only predicts satisfactorily their effective charge, but 
also reproduces accurately the anisotropy of their po- 
tential. The latter property, embodied in the function 
f(na, 9) is a key feature of screened electrostatic interac- 
tions and may have non negligible -and hitherto largely 
unexplored- consequences. It is responsible for the rich 
phase behaviour and orientational ordering of colloidal 
molecular crystals |29| . Its effects on the phase behaviour 
of clays, especially at moderately to high salt concentra- 
tions where large energy barriers f(8 = it/ 2) — f(9 = 0) 
are observed, will be the subject of future work. 

The authors acknowledge fruitful discussions with J.J. 
Weis, B. Jancovici, F. van Wijland, M. Aubouy and H. 
Lckkerkerker. 



APPENDIX A 

Mixed boundary value problems are rather frequent in 
electrostatics but also in diffusion and elasticity prob- 
lems [3(], El (conduction of heat, diffusion of thermal 
neutrons, punch or crack problems etc.). They may be 
encountered in hydrodynamics as well |32j |. They gen- 
erally arise whenever a potential is prescribed over part 



of a boundary whereas its normal derivative is specified 
over the complementary part. The theory of dual inte- 
gral equations turns out to be a powerful tool for such 
situations. In this appendix, we give more details about 
the problem of finding the solution A{k) to equations @ 
and (jSJl. To begin with, it is convenient to recast them 
in the form 



Iq Ju 2 + (no,) 2 
g(u) Jo(xu) du 



Jo (xu) du 



1 for x < 1 



for x > 1 



(Al) 

where dimensionless quantities have been introduced: 
u = k a, x = p/a and g(u) = y/v? + (na) 2 A(u)/$ a . So- 
lutions of the previous equations for k = (no salt case) 
have been derived by Titchmarsh [3j|. The procedure, 
based on rephrasing the dual integral equations by means 
of some invertible linear operators gives the solution 



A K=0 (u) = 



2$o sintt 



(A2) 



and the corresponding potential is that which leads to 
equation l|14fl for the charge. A generalization of Titch- 
marsh's method has been proposed by Sneddon [34| for 
dual integral equations of the type 



u 2a (1 + ^(u)) g(u) J v {xu) du = 1 for x < 1 



g(u)J v (xu) du 



where ui is an arbitrary function. Equations i|A3|) and 
(|A1|) can be made equivalent by taking 



= for x > 1 
(A3) 



uj(u) 



^Ju 2 + (na) 2 



(A4) 



with a — 1/2 and v = (for uj(u) — 0, we recover the no 
salt case). The problem at hand -that fits into the gen- 
eral framework of [34[- may be reduced to that of solving 
a Fredholm equation of the second kind. Following Sned- 
don, we write equations i|A3|l in the form: 



(A5) 



S_i 1 [(l + uj(u))g(u)/u, x] — 2/x for x < 1 



So,o[g(u)/u,x] =0 for x > 1 

where S Q .^ is the modified Hankel operator defined by: 
S a</ 3[\(u),x] = S a ,f3\(x) 



u 13 X(u) 32a+p(xu) du 



We then introduce the function h through 

g(u) = uS ih(u). (A6) 

This function h(u) is the central object in the present 
procedure. After some cumbersome algebra based on the 
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pro pert ies of the modified Hankel operators (for details, 
see |34|); we find that the function h 2 = u h(u) defined on 
[l;oo] vanishes while its counterpart hi defined on [0; 1] 
is the solution of the following Fredholm equation of the 
second kind: 

hi{x)+ ( hi{u)K{x,u)du = (A7) 



The kernel K (x, u) is defined in the whole xu-plane by: 



K(x, u) 



Ka(Ii(Ka\x — u\) — Li(kcl(x + it))) (A8) 



and I\ and L\ denote respectively the modified Bessel 
function and Struve function of the first kind, of order 
one. The weight function g(u) appearing in equation 
(|A1|) is recovered by means of the integral: 



g(u) = yju 2 + {no) 2 A(u)/$ 



cos (ut) hi(t) dt. 



(A9) 
(A10) 



The potential finally follows from equation The Fred- 
holm equation (|A7(I is solved numerically by an iterative 
procedure, starting by an constant initial guess for h\. 



Ka=0.1 
Ka=l 
Xa=2 
Ka=5 



0,2 



0,4 0,6 
x 



FIG. 18: Plots of hi(x) [solution of the integral equation IA7I 1 
versus rescaled distance x for different values of Ka. For Ka — 
0.1, hi is very close to its no salt limit 2/ v / 7r — 1.13. 



In the limit K = 0, K(x, u) = so that, from Eq. I|A7() . 
h\ = 2/y / 7t. The function g(u) in (|A10|) follows immedi- 
ately: g(u) = (2/-7r) sin it, which is fully consistent with 
(|A2|) . Finally note that the functions hi(x), plotted in 
Fig. |^| for different values of Ka, are related to the sur- 
face charge density of the disk a{x) through the equation: 



a{x) 



h x (x) 



1 



(All) 



The function h\ can be well approximated by a quadratic 
polynomial. The corresponding charge cr(x) may then 



be used to compute analytically approximated effective 
charges, anisotropy functions and weigh functions A(u). 
This procedure is an alternative to the ansatz proposed 
in equation l|15l) . the latter being more suited for an an- 
alytical treatment. We finally note that from Eq. (|A11|) 
and the regular behaviour of hi observed on Fig. 1181 for 
x — > 1, the surface charge diverges near the edge of the 
disc like (1 — x)~ x / 2 , see appendix iBl 



APPENDIX B 

We show here that the LPB surface charge distribu- 
tion <j{p) -arising from the condition of constant surface 
potential on the disc- exhibits in an electrolyte the same 
edge effect as in vacuum, where it diverges as (a — p) -1 / 2 
when p — > a~ [2(j . 



I— 



B = ti/2 B = 7i 



B = 3tc/2 B = 2tc 



FIG. 19: The wedge geometry. The conductor is represented 
as the shadowed part and its boundaries by the thick line. 

We consider the more general problem of a wedged- 
shape conductor with angle f3 (see Fig. IT^|) . The imposed 
potential is denoted cf>o. Being interested in the behaviour 
near the sharp edge where the situation is of cylindrical 
symmetry, we introduce the cylindrical coordinates (r, 
xp) in a plane perpendicular to the apex of the wedge. 
With rescaled distance r = nr, we look for solutions of 
LPB equation (01 



r r — — 

dr \ dr 



cV 2 



in a form with separated variables 

<j>(r,ip) = R(r)^(ip). 



(Bl) 



(B2) 



With the boundary condition 4>(r, ip = 0) = </>(r, ip — 
(3) = <fio for all r in the vicinity of the wedge (i.e. r close 
to 0), the solution reads: 



R{7)=A v I v {r) 
<&(V>) = ol v sin(z/'i/')- 



(B3) 
(B4) 



Here, v = rnr/P where n G N; A v and a v are arbitrary 
constants. In the vicinity of the wedge, the potential 
therefore takes the form 



(j)(r,ifi) = (fi 



oo 

E 

n=l 



A n sin 



mrip 



In*/(i<r)- (B5) 



The dominant term when r — > corresponds to n = 
1, and scales like r'' 15 . The associated surface charge 
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behaves like a(r) oc r 71 ^ 1 , hence like r x l 2 near the 
edge of a disc ((3 = 2tt). This result has been used to 
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